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We study the spacetime-reduced (Eguchi-Kawai) version of large-iV QCD with nonzero 
chemical potential. We explore a method to suppress the sign fluctuations of the Dirac deter- 
minant in the hadronic phase; the method employs a re-summation of gauge configurations 
that are related to each other by center transformations. We numerically test this method in 
two dimensions, and find that it successfully solves the silver-blaze problem. We analyze the 
system further, and measure its free energy F, the average phase 9 of its Dirac determinant, 
and its chiral condensate (^if>). We show that F and (ipip) are independent of /i in the hadronic 
phase but that, as chiral perturbation theory predicts, the quenched chiral condensate drops 
from its /i = value when \x ~ (pion mass)/2. Finally, we find that the distribution of 9 
qualitatively agrees with further, more recent, predictions from chiral perturbation theory. 
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I. INTRODUCTION 



It is of great interest to calculate properties of four dimensional QCD at low temperatures 
and large, but not asymptotic, chemical potentials. Due to the 'sign problem' that inflicts 
euclidean lattice Monte-Carlo simulations such calculations seem currently unrealistic [l| (for 
recent progress in this see Ref. j^j]). Specifically, the fluctuations in the phase 6 of the Dirac 
operator's determinant, which become particularly strong once the chemical potential /x grows 
beyond half the pion mass m n /2, cause the simulations to fail. At a first glance this seems 
surprising because, at least for low temperatures T <C Aq CD , one expects physics to depend 
on fj, only once fi > wib/3, with mg the lightest baryon mass. Thus, for m^/2 < \x < tub/S, 
physical observables are approximately independent of \i but 9 is strongly oscillating. The 



apparent tension between these two facts was termed the 'silver-blaze' problem and for 
small /i it can be studied with chiral perturbation theory [4]. The sign and silver blaze 
oroblems were also recently studied in scalar field theory using complex Langevin dynamics 
5| and in fermionic theories with four-fermi interactions [6|. 

In this paper we analyze various aspects of the QCD sign problem. In particular, we ask 
how the sign problem behaves once we project the Dirac determinant to be neutral with 
respect to the Zm center symmetry of the SU(N) group. Our current paper is not the first 
to discuss this 'Z^-averaging' — see for example Refs. [g] - but we are not aware of studies 
that explicitly checked its effect on the sign fluctuations of the determinant in the confined 
phase of QCD. 

Our modest computational resources lead us to approach the problem in its large-iV 't 
Hooft limit and to utilize the space-time reduction of the planar theory. This reduction in 
degrees of freedom was discovered by Eguchi and Kawai in the seminal Ref. [q| whose content 
is the following: under certain conditions, infinite-volume lattice large- N QCD is equivalent, 
on all distance scales, and in certain sectors of its spectrum, to its 'volume-reduced' version. 
The latter is nothing but lattice QCD defined on a single lattice site. 

The simplest implementation of the Eguchi-Kawai (EK) equivalence works only in two 



euclidean dimensions 
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prescription, like those of Ref. 12 



e in four dimensions, some extensions of the original EK 



13], are expected to be successful. The numerical work in- 
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volved in using these extended prescriptions is comparably demanding, and their applicability 
is still being tested numerically [ljl [lsj . These facts lead us to focus on the two-dimensional 



system in this paper. Provided that the constructions of Refs. |12l . [13] survive their numerical 
tests, we do not see any obvious numerical or conceptual difficulties in extending the present 
approach to four dimensions. 

Below is the outline of our paper. We begin in Section [TT] by reminding the reader what 
are the validity conditions of the EK equivalence and discuss what they imply for our nu- 
merical calculations. We then present the construction of the two-dimensional EK theory in 
Section and discuss how we simulate it and which observables we measure. We formulate 
Zjv-averaging in Section IIVI Then, in Section [V] we discuss the connection of the sign and 
silver-blaze problems to the signal to noise ratio (SNR) of physical observables calculated with 
a Monte-Carlo simulation, and the way averaging can potentially solve these problems. In 
Section IVII we present our numerical study of the sign problem and show the average sign of 
the determinant as a function of /i. The tests of Z^- averaging are presented in Section IVIII 
We then show results from measurements of various physical observables in Section IVHH and 
in Section HXl we present the distributions of 9 as a function of /i, and the way the eigenvalues 
of the Dirac operator scatter in the complex plane for different values of \x. In Section |X] we 
summarize our study. 



II. EGUCHI-KAWAI SPACE-TIME REDUCTION AT NONZERO DENSITY 

There are two issues one should be aware of when one studies nonzero density with space- 
time reduction. We discuss these issues in Sections III Al and III Bl where we show that the 
method of EK space-time reduction can be used only within the hadronic phase. Despite that, 
this method is useful to understand the silver blaze and sign problems of QCD at nonzero fi, 
and we emphasize this point in Section III C[ 
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A. The importance of translation symmetry 



In Ref . 



161 ] we studied the space-time reduction of large- N QCD in the presence of baryons. 



To have full analytic control we focused there on the two-dimensional version of QCD - the 't 
Hooft model. We found results which are consistent with the general validity conditions that 
non-perturbative large- N equivalences require. Specifically, we saw that the main validity 
condition for volume reduction (or the 'Eguchi-Kawai equivalence') in our context is that 
translation symmetry is intact in the infinite- volume field theory. 1 To understand this recall 
that the EK theory is the single-site version of QCD. Thus, it is a projection of QCD that 
removes degrees of freedom with nonzero momentum. As such, it is clear that using the EK 
equivalence does not make sense when the vacuum of QCD spontaneously breaks translation 
invariance and is characterized by condensates which carry nonzero momentum. The question 
of whether translations are a symmetry of the QCD ground state at nonzero \i is a dynamical 
one, and in Ref. JlT] we showed that in two dimensions, for any quark mass, and for a single 
flavor, this symmetry breaks spontaneously for \i > ttib/N (for earlier results concerning only 
the chiral limit see Ref. 18|). Simple arguments lead one to expect a similar phenomenon 
in higher dimensions (see discussion in section IV of Ref. This means that the single 

site theory that we study in the current paper is equivalent to the infinite- volume field theory 
throughout the hadronic phase which, in two dimensions and in the single flavor case, is 
bounded by fi < m B /N. Beyond that point the theory we study here can be thought of as a 
complicated matrix model. 



B. Implications of lattice saturation 

As any lattice field theory, volume-reduced QCD can be expected to be influenced by lattice 
artifacts at sufficiently large values of //. In particular, for lattice spacing a and chemical 
potential fi ~ 0(l/a), the density is at the cutoff scale, and the lattice is saturated with 
baryons - there is an 0(1) number of baryons per site. This regime, in which the Pauli 

1 Another condition is that the Zn center symmetry of the theory is intact, but this is dictated by the leading 
gluon dynamics and not by the fermions. 
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exclusion on each site is saturated or nearly saturated, is usually referred to as the lattice 
saturation regime and is governed by lattice artifacts. 

What does the saturation phenomenon imply on a single-site theory? In d + 1 dimensions 
the spatial volume V of the single-site theory is a d , and as long as /x is smaller than a critical 
value /i c , there are no baryons in the vacuum (as mentioned above, for d — 1 and with a single 
flavor, baryons repel and \i c = rns/N 17|, [l9j). But when \i > fi c the hadronic phase makes 
way to a phase that accommodates baryons, and the baryon number B is then at least 1. 
Importantly, this makes the baryon density, B/V, of 0(l/a d ) and so at the cutoff scale. 

Therefore, we see that we cannot accommodate physical nonzero densities on a single site: 
when the baryon number is nonzero, it is at least one, the density of the single-site theory is 
at the cutoff scale. 



C. Usefulness of space-time reduction at nonzero fj,. 

At this point an obvious question comes to mind — what is the usefulness of the EK 
single-site theory at fi > and why do we wish to study it in this paper? The answer is that 
we wish to understand the sign and silver blaze problem of the hadronic phase. In that phase 
fi is nonzero, but the density is zero and so we can use the EK theory to analyze the theory 
there. 

Also, despite the fact that in the saturation regime (outside the hadronic phase), the EK 
theory is not equivalent to the field theory, we can still ask how the EK model behaves 
there, and see whether we can learn something about the general properties of the saturation 
regime (as we alluded to above, this regime appears also in the standard, infinite-volume, 
'non-reduced' lattice field theory). 

III. LATTICE DETAILS OF THE EGUCHI-KAWAI THEORY AT NONZERO fj,: 
THE ACTION AND THE SIMULATION ALGORITHM 

The EK theory or 'volume-reduced QCD' is a lattice gauge theory defined on a single 
lattice site. Since our numerical efforts are focused on two dimensions, we define below the 
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construction of the EK theory only in this case. 



A. Definition of path integral 

The path integral is 

Z = J DU J DijDij exp (S YM + S F ) , (3.1) 

S YM = 2Nb Re Tr (u x U 2 U\ U ] 2 ) , (3.2) 

SWac = 4>Dip, (3.3) 

D = m + i 7l (u x e A - e-* U\) + \l2 (u 2 - U ] 2 ) , (3.4) 

Here U\^ are SU(N) matrices and DU = dUydU^ with dlli^, the Haar measure on SU(N). 
ip is a two dimensional Dirac spinor of the 'naive' fermion type and it transforms in the 
fundamental representation of the gauge group. Thus, it corresponds to 4 x Nf Dirac fermions 
in the continuum. The 2NfN x 2NfN matrix D is the euclidean Dirac operator of a lattice 
gauge theory on a single site. The quark mass m and chemical potentials fi are related to 
their dimensionless lattice quantities rh and p, in the usual way by the lattice spacing a. In 
two dimensions the gauge coupling g has dimensions of mass and the standard dimensionless 
lattice coupling (3 is defined by 

2N , . 

P = ~T2- 3.5 

In the 't Hooft limit g 2 scales like 0(iV _1 ) and so it is useful to define a new lattice coupling 
b = P/ (2N 2 ) that scales like 0(N°) in the large- N limit - this is the coupling that appears in 
Eq. (13.21) . In terms of the dimensionful 't Hooft coupling A = g 2 N, the coupling b is given by 

1 



a 2 X 



(3.6) 



All dimensionful quantities like the quark mass m and the baryon chemical potential /i will 
be fixed in units of A when we take the continuum limit b — > oo. Also, to make a connection 



with Ref. 20], we define the parameter 7 that controls the quark mass 



2 

m 



7 = Tr— • (3-7) 
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In terms of b and 7 the lattice quantities are given by 

m 1 FY~ 

m = am= 7s7rvs- < 3 - 8 > 

^ v = 7xTb (3J) 

It is useful to note the following results on the spectrum of the 't Hooft model in the large- 
N limit: The pion mass at zero chemical potential for small and large values of 7 is given by 
(see Ref. [20 ] ) 

\ 1.08 7 1/4 + O( 7 1 / 2 ) 7«1, ,_ ln , 
— = ~ < (3.10) 

V\ [ 1.13 7 1/2 + 0( 7 - 2/3 ) 7>1, 
which is expected to hold for any number of flavors Nf that obeys Nf <C N. Next, the lightest 
excitation with a nonzero baryon number in the single-flavor case is the baryon whose mass 
in the chiral limit is 19 1 



^ = lx^. (3.11) 
N n 2 K J 

Indeed, in two dimensions, both m w and m# vanish for massless quarks! For a non-zero quark 
mass, however, ms/N > m n /2 and we have the intermediate range, m n /2 ;$ jj, <, uib/N, 
where the theory has the same properties as those at \i = 0, but where the sign fluctuations 
are strong — just like in the case of physical four dimensional QCD. 



The most recent calculation o 



t 



iems/iVinl + 1 dimensions 



'or Nf = 1), away from 
the chiral limit, was done in Ref. 17[. The previous study in Ref. [19[ had certain ingredients 
missing, whose effect was calculated in Ref. and was found to be small. In any event, 
using the methods of Ref. [17| , we evaluated the baryon mass for the quark mass studied in 



the current paper and we find 2 

m B 



0.758 at 7 = 1. (3.12) 



nVx 

We note that both the estimates for in Eqs. ( 13 . 1 2 j) and (13. lip are for Nf = 1 and we are 
not aware of any extensions to general Nf. Since we define the EK theory with naive fermions, 
for sufficiently small lattice spacings we actually have four physical Dirac fermions. In the 



2 This estimate was obtained in the Hamiltonian formalism at a spatial lattice spacing of a^/\ ~ 0.31. The 
lattice spacing corrections are expected to be small for this value of a [l7|. 
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absence of any estimates for the baryon mass in the four-flavor case, we proceed by assuming 
that the Nf dependence of mg is weak. 



B. Calculation strategy 

Expectation values are measured with Monte-Carlo simulations that use the Yang-Mills 
action to generate gauge configurations. This means that (O) is approximated by the sum 
over gauge configurations that represent the YM ensemble: 

1 M 

(O) * s g Oc (3.13) 

Here M is the number of configurations. Our measurements can be therefore referred to as 
performed via 'reweighting' because we generate the gauge configurations with the YM action, 
and reweight in the fermion determinant. 

This reweighing would not be necessary had we focused on \x = since the large- N theory is 
naturally quenched there (at least at nonzero quark mass); in practice this would be reflected 
in a Monte-Carlo by the fact that the 'fermionic force' in a Hybrid Monte-Carlo is 0(1/ N) 
suppressed compared to the gluonic force. At \i > 0, however, reweighting in the phase of the 
determinant is crucial (as our results would show). 

By reweighting from the YM ensemble, we use the fact that the important gluonic config- 
urations at large- N are chosen by the YM action, and that the fermions do not back-react 
to the gluons. This absence of back-reaction is characteristic of the large- N limit and can 
be justified in various approaches (like Hamiltonian coherent states) but we will not discuss 



this issue here (see instead the forthcoming Ref . [30[ ) . We should note that the lack of back- 
reaction at large-iV persists at nonzero /i, but that this does not mean that the quenched 
approximation should work there (despite the common lore, the quenched approximation is 
not equivalent to removing the back-reaction of quarks on gluons from QCD - more details 



will be given in Ref. [301] ) . 

In the context of standard lattice QCD studies, the reweighting procedure is known to be 
highly susceptible to significant systematic errors when fi and the lattice volume are large 
enough [l] . These errors can be diminished, however, if one scales the number of sampled field 



configurations exponentially in the size of the matrix D. In four dimensional 3-color QCD one 
has dim(D) = 12 x (lattice volume), and so going to the thermodynamic limit is not realistic 
with current computational power. Similarly, in our case we have dim(D) = 2 x N and our 
computational resources allow us to study only N < 60. 3 



To simulate SVm we used the heat-bath algorithm introduced in 211 ] . Measurements of 
fermionic quantities are typically separated by 100 — 1000 full updates of the model, 4 and thus 
are expected to be uncorrelated. To check this we estimate errors with a jackknife procedure. 
Below we list the fermionic quantities that we measure for most values of b, N and \i. 

a. The fermionic contribution to the free energy, F. This was obtained from 

F = 1 log Z QC D = ^ log (det D(jjl)) . (3.14) 

b. The average sign of the determinant given by 

<-w>-M£a))- (3 - 15) 

c. The quenched quark condensate which we define to be 

quenched 

(tr (zrV)))- (3.16) 

d. The physical unquenched quark condensate 

/7,\ (tr (D" 1 ^)) xdetD) 

(M) unquenched = — (3-17) 

e. The distribution of the angles 9 and a defined via 

det£> = | det D | e ie , (3.18) 
^ = ^ e ia . (3.19) 



3 Our errors modestly increase with N, and this is because we do not actually increase our statistical sample 
exponentially. Nevertheless, we find that N = 40 and 60 are already close enough to the large- N limit and 
in most cases have errors that are sufficiently small for our purposes. Indeed, this is how large- N reduction 
plays an important role in our study. 

4 By 'full update' we mean updating all the N(N — l)/2 SU(2) subgroups of each of the two matrices Ux s- 



S 



f. The way the eigenvalues of D scatter in the complex plane. 



Finally, for fi = and for each value of b and N we also measured the pion propagator of 
momentum q along the "2" direction using the Gross-Kitazawa momentum injection method 



22| | (that was also used in Ref. 23]). This means that we measure 

G w (q) = (tr (L>- 1 (C/ 2e ^ 2 ) T5J D- 1 (C/ 2 e-^ 2 )75)) , (3.20) 



and extract the pion mass by measuring the effective mass m e s defined by 

r 3 [log (Refdqe^G^q))] 
m eff = - hm . (3.21) 

We calculated these observables for chemical potentials in the range fi/Vx e [0,3]. 
Throughout our calculation we fixed 7 = 1 and so we expect m n /\/\ ~ 1.2 from Eq^_ (I3.10p 
and m^/y/X ~ 1.5 from the continuum extrapolations of the numerical results in [20(. We 
perform lattice simulations at b = 0.6, 6.0, 10.0. The number of field configurations that we 
used for each choice of N and b is given in Table E We also list the pion mass in each case 
as measured from Eq. (I3.2ip . For brevity, we set ^m n /y/\ = 0.65to;o!5 to encompass all the 
values of m 7r /2 from Table [B We present the results of our study in the next sections. 

The numerical routines we used were defined with double precision, and to check whether 
our calculations are sensitive to this (especially our evaluations of the determinants) we recal- 
culated the free energies and average signs with quadruple precisions and confirmed that our 
results remained the same. 



IV. DEFINITION OF Z N - AVERAGING 

In this section we define Z^- averaging and briefly explain its logic. This averaging relies 
on the Zn symmetry of the YM part of the action and of the Haar measure, and it is defined 
by the following prescription. 

For each gauge configuration C, characterized by the gauge fields (Ui,U 2 ), perform the 
replacement 

1 N 

O e (U u U 2 )^(O c ) z = — Oc^e^^^e^/"), (4.1) 

fel,fc2=l 

9 



N 


b 


no. of configurations 


m n /V\ 




0.6 


16000 


< 1.20 


10 


6.0 


20000 


1.24(12) 




10.0 


20000 


1.24(13) 




0.6 


12000 


< 1.24 


20 


6.0 


120000 


1.25(5) 




10.0 


120000 


1.33(6) 




0.6 


1000 


< 1.29 


40 


6.0 


105000 


< 1.23 




10.0 


101000 


< 1.40 


60 


10.0 


50000 





TABLE I: Details of runs used to map the phase diagram in /z. The number of full model updates 
was of O(10 6 — 10 7 ) depending on the value of N. Not all configurations were used to estimate the 
meson masses. The cases in which we only give an upper bound on are those in which we did not 
observe the asymptotic mass plateau in the pion propagator. The reason for this is either that m, 
was too large in lattice units, and so falls into the noise quickly (we saw this for b = 0.60), or that 
one needs to calculate the pion propagator at separations x which are too large for our resources to 
accommodate. 

and only then average over the field configurations U\$. We distinguish between the results 
obtained with and without this Z^- averaging by adding the subscript Z to the former, when 
necessary. In four dimensions Eq. (14. ip is generalized to have four sums over ^2,3,4 £ [1 , A/"] . 

Let us explain the logic behind this proposal. We expect that the partition function will 
depend on /1 in a very special way; for a U(N) group there should be no //-dependence since 
there are no gauge invariant charges that couple to \l. For an SU(N) group, however, there 
can be /i-dependence, but only through the combination /xN since the only gauge invariant 
excitations that couple to fi are baryons, and these have charges that are integer multiples of 
N. 

How can we anticipate these results from the path integral in Eq. ( 13. ip ? Expand detD in 
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Eq. ( 13.1 p in the worldline approach [7f . It is then easy to see that the //-dependence we discuss 
above arises because both the measure of the path integral and the YM action are symmetric 
under the center of the gauge group. Specifically, the \i dependence of the determinant comes 
in through terms of the form 



e 9 " x [] (trUt) P \ and e 



x 



(4.2) 



with q = J^ikiPi — HiKv'ii but the center symmetry allows q = for U(N) and q/N = 
integer for SU(N)\ Z^r- averaging enforces the center symmetry on each gauge configuration, 
and automatically makes the path integral depend on fi as prescribed above. In fact, since it is 
only U\ that appears in Eq. (14. 2p . we can also try and do a partial averaging by re-summing 
gauge configurations that differ by a center transformation only in the time direction. 

N 

X 



1 N 

O c (U u U 2 ) -+ (O c ) z _^ rtml = ^Y.°c (t/i e^l\ U 2 ) . (4.3) 



In our two-dimensional case this saves a factor of N in the measurements of observables (and 
jyo-i f or d space-time dimensions). Numerically we see that, at least in the hadronic phase, 
such a partial averaging is almost as good as the full Zn averaging (see Section IVlII) . and 
so we use it when our resources are insufficient to perform the latter (for example when we 
calculate the unquenched quark condensates for SU(4Q) or the average sign for SU (60)). 



V. CAN Z N AVERAGING SOLVE THE SILVER-BLAZE AND SIGN PROBLEMS? 

As explained above, one reason that we expect Z^ averaging to be useful is that it makes 
the anticipated dependence of the partition function -Zqcd on A* manifest on a configuration- 
by- configuration basis — it is a noise reduction technique. To make this expectation more 
precise we need to discuss the source of the statistical noise in the calculation of -Zqcd for 
/i > (and its relation to the sign problem). 5 We do so in Section TV Al and then move on to 
Section IV Bl where we discuss how such noise can be reduced when one uses Z^ averaging. 



For a previous useful discussion on see Ref. [24 1 
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A. The sign problem, the signal to noise ratio, and pion physics. 



The sign problem and the signal to noise ration (SNR) of observables in QCD are closely 
related. To show this, let us discuss the SNR of observables in a Monte-Carlo that generates 
gauge configurations for the 'Phase-Quenched' (PQ) ensemble. The measure of this ensemble 
is proportional to the absolute value of the Dirac determinant. Specifically, if we write 

detD = \detD\e i9 } (5.1) 

then, when calculating the expectation value of an observable O with gauge configurations of 
the PQ ensemble, one needs to treat the phase e %e as part of the observable. This means that 
the QCD average of O becomes 

<^<gH (5.2) 

It is useful to note that we can define the PQ average (, )pq through the Yang-Mills average 
(, )ym 

( ° )PQ = (|detD|) YM - (5 ' 3) 
From Eq. (I5.2p we see that the SNR of any observable gets a contribution from the SNR of 
(e i6, )pQ = (cos6')pq (here we assume charge conjugation). The SNR of (cos6 i )pq is given by 



SNH2, S ^9)!. (5.4) 
(cos 2 #)pq 



In terms of (, )ym, Eq. ( 15. 4 p becomes 



QNR PQ = ((cos0|det£>|) Y M) 2 = ((det£>) YM j 

cose (cos 2 9\ det _D|)ym(| det _D|)ym (cos 2 9\ det D\) YM (\ det D|) YM ' 1 ' ' 

and it is now clear that if there is no sign problem, and cos# ~ 1 for the majority of field 
configurations, then SNR CO g ~ 1. 

It is instructive to show how pion physics tends to appear when there is a sign problem. 
This can be seen by assuming that a sign problem makes cos 2 Q —\ and turns Eq. ( 15. 5D into 



\2 
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Now, for simplicity, consider the case of two degenerate flavors and note that 



( 



(detDf>YM 



) 



2 




6 



~ 2 



(IdetDiHvM 



2e +2V/T(f I -f B ) 



(5.7) 



Here D 1 is the Dirac operator of a single flavor and fs and // are the free energy densities 
of two-flavor QCD in the presence of baryon number and isospin chemical potentials that are 
both equal to /i. To proceed, note that baryons do not contribute to /# for (j, < mg/JV (we 
always restrict to low temperatures) which gives 



In contrast, // changes when \x > m n /2 (the iso-spin system goes through pion conden- 
sation above this value of //). These facts, together with // < which is true because 



B. Two scenarios for the numerical cause of the sign problem and how Zjy averaging 
can solve it 



dated sign problem from the point of view of the worldline approach. In particular, these 
problems reflect sign fluctuations in det D that are caused by the contribution to (det D) from 
Polyakov loops that wrap around the temporal direction. From this point of view, there are 
two scenarios for the cause of the sign problem: 

1. Scenario no. 1: The sign problem is caused by fluctuations of Polyakov loops that wrap the 
torus k times with k/N / integer. 

Polyakov loops whose winding number around the temporal direction k is not an integer 
multiple of 3 (or iV in SU(N)) are unphysical; by the center symmetry these nonzero iV-ality 
worldlines have zero expectation value and so they do not contribute to the numerator of 



/b(/x) = /b(0) for fi<m B /N. 



(5.8) 




We can understand the problem of the exponential suppression of SNR^ e and the asso- 
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Eq. ( 15. 7p . This is despite the fact that they appear in det D\ with enhancing fugacity factors 
of e k ^ T . 

While the unphysical worldlines do not contribute to (det Df) YM , they do contribute to 
its noise (|det i^i | 2 ) appearing in the denominator of Eq. (15. 7ft ; despite the fact that their 
magnitude can be small, their fugacity factors e k>J, ^ T can be large, and they can cause significant 
fluctuations in the overall fluctuating phase of the determinant. The magnitude squared of a 
fc-worldline contribution to det D\ is a worldline configuration with k worldline-antiworldline 
pairs which can be interpreted as a worldline of k pions (this can be seen explicitly from the 
worldline expansion of the |det | 2 ). Therefore, this magnitude squared can be estimated 
by e~ fcm " r//T , and we can now anticipate that det will have significant noise when e fcM//T x 
e -k m7T /2T ^ i 0T w hen jj > mjr/2, as we saw in Eq. (15. 91) . 6 

If the scenario we describe above is the one causing the sign problem then it is an optimistic 
scenario: worldlines that wrap the euclidean time torus k times with k/N ^ integer can be 
removed from the noise on a configuration-by-configuration basis using Zn averaging. This 
would solve the silver-blaze problem since it would remove the phase fluctuations of the Dirac 
determinant in the hadronic phase. 

Importantly, however, our focus on the unphysical worldlines in the current scenario 

assumes that the physical, 'baryonic', worldlines with k/N = integer, are suppressed on 

a configuration-by-configuration basis (assuming, for example, that they are weighted by 
e -km B /Ty If that ig the 

case then the only sign fluctuations that would survive aver- 
aging would occur when the baryonic worldlines contribute to (detDf), i.e. when /x > tub/N. 
In the next scenario below we ask what happens if this is not the case. 

2. Scenario no. 2: The sign problem is caused by Polyakov loops that also have k/N = integer. 

Unfortunately, the physical, 'baryonic', worldlines, for which k is an integer multiple of N 
(and therefore are insensitive to Zn averaging), can have a sizable contribution to the sign and 

6 In principle we should replace km^ by the free energy of a system with isospin number equal to k, but the 
binding energy of k pions is suppressed at large- AT. 
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SNR problems even at small fi. Naively one might think that these worldlines are suppressed 
so long as fi ^ uib/N (see end of previous section) but this is not necessarily correct. To see 
why, let us focus on the single-baryon contribution to (det-D), and denote it by (det D) B=1 . 
The latter reflects the contributions of all worldlines that wrap the euclidean circle N times 
and so it comes with a fugacity pre-factor of e N ^ T . Our main point here is that there is no a 
priori reason why these contributions should be suppressed by e _ms//T on a configuration- by- 
configuration basis. Quite the contrary, it is possible that there are some gauge configurations 
where these contribution are of O (p- Nm */ 2 \ — here we view the weight of a single quark line 
as the weight of half of a pion. 7 

Thus, we propose that (det D) B= ^ evaluated on a single gauge configuration would have 
the following generic form reflecting the fact that its noise (or absolute value squared) has a 
combined origin of both a baryon-antibaryon pair or iV pions. 

(det D) B=1 « e^»-™ B )IT + c x ^/a)/r (5 10) 

Here the (generally complex) number c depends on the dynamical details of the theory and on 
a configuration- by-configuration basis can be of an 0(1) magnitude but strongly fluctuating. 
In particular, the fluctuations in c can average the second term in Eq. (I5.10p to zero, such 
that the gauge-configuration average will make ((detD) B=1 ) (that determines the baryon 
mass) associated only with the iV world-lines that bind into a baryon. 

From Eq. (15.10!) we see that if we send T — > 0, then the second term is the dominant 
one (since m n is always smaller than 2ms /N). This means that when fi > pions will 

proliferate in the noise of detD and there will be a severe sign problem. In contrast, if we 
keep T finite (but low), and if |c| happens to be numerically very small, or more precisely if 

\ C \ < e W2m,- mB )/T^ ( 5 _ n) 

then the pionic term can become irrelevant. In that case (det D) B=1 will be of the order 
of e~ ms//T for most of the gauge configurations, and all the fi dependence of det D will be 

7 Indeed, the noise in the evaluation of (det D) B=1 is given by the magnitude squared of (detD) s=1 , and 
can describe either a baryon-antibaryon pair or an iV-pion system (both can be described by N worldline- 
antiworldlines pairs). 
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exponentially suppressed if fx < m B /N. Then detD will behave like it does at fx = 0, 
i.e. there will be no sign problem. 

It is Eq. (15. lip that determines whether averaging solves the sign problem or not. 
Similar arguments lead the search for optimized baryonic wave functions that couple minimally 
to pion physics, thus improving the signal to noise ratio of baryon correlation functions 26] . In 
that context, one searches for a 'golden window' in the euclidean time separation of baryonic 
correlators, where their SNR is not exponentially small. Here, Eq. (15.111) defines a golden 
window in the temperature T. 

In the following sections we present results from numerical simulations, some of which 
done in order to determine whether Z^ averaging indeed removes the sign fluctuations in the 
hadronic phase and solves the associated silver-blaze problem. Put differently we aim to see 
whether it is scenario no. 1 or no. 2 that takes place in our system. 



VI. RESULTS : THE SIGN PROBLEM 



In this section we present the numerical analysis we performed for the averaged sign. We 
begin by presenting the way the average sign of det D behaves as a function of fx for SU(40) 
and different values of b in Fig. [IJ The vertical (gray) solid band in the figure denotes our 
estimate for m 7r /2 0.65a/A (see above). The dashed vertical (red) line denotes the baryon 
mass ffig (divided by N) for the single-flavor case. It is higher than m n /2 by only about 
16% (this proximity is special to 1 + 1 dimensions). As we emphasize in Section [Til the EK 
theory that we study is equivalent to the 2d gauge theory on an infinite volume only below 
uib/N. We nevertheless present here (and below) the results for larger values of fx, since for 
these values of fx the model has a saturation behavior which is accompanied by an interesting 
behavior in the average sign. Specifically, from Eq. (13 .4p we see that when the EK theory is 
saturated with baryons, it has no sign fluctuations because the term governs the behavior 
of D and one has 

L> ~ C/ie^ (l + 0(e~ A )) . (6.1) 
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Average sign, SU(40) 




FIG. 1: The average sign, (cos(0)) for SU(4Q) and b = 0.60,6.00,10.00. 
Since U\ G SU(N) this means that 

det£) saturation ^ ( g ^ 

and that (cos8) — > 1. Looking at Fig. [Q we see that this is happening at /isaturation/'s/A ~ 1.5 
for b = 0.60, and starting to happen at /i sa turation/"\/A ~ 3.0 for b = 6.00. The increase with 
b of /isaturation is expected since in the continuum limit, b — > oo, the saturation goes away 
completely (recall that this saturation reflects the fact that the baryon density on a lattice is 
bounded from above by the 'saturation density', and the latter is of 0(l/a d )). The suppression 



in the sign problem in the saturation regime was also seen in the work of Ref. 



27] 



Another interesting fact is that (cos 6) drops from 1 at /i ~ rn n /2 - as chiral perturbation 



theory predicts (see, for example, Ref. [4] and its references). Importantly, we point out that 
above fi/VX ~ 1.0 the average sign is already very close to zero and so there is a severe sign 
problem there. 

We now turn to study the way (cos 8) changes with N. Since the large iV limit is the 

17 



thermodynamic limit in the EK theory, we expect the sign problem to become worse as N 
increases. This expectation is confirmed in Fig. |2]where we plot the average sign for b = 10.00 
and for the gauge groups SU(10), SU (20), SU(40) and SU (60). We see that while there is no 

Average sign, b=10.00 




FIG. 2: Comparing the average sign of SU(N) with N = 10,20,40 and 60 for b = 10.00. In a 
dash-dot (blue) line is the average sign of the epsilon regime. 



serious sign problem in SU(10) for any value of fi, for SU (40) the average sign approaches zero 
at around fi/y/X — 1.25, and for SU (60) the average sign is zero already at fi/y/X — 0.75. The 
way the average sign decreases with N is expected to be exponential if \x > m n /2. This can 
be easily argued if a different definition of an average sign is used (for example, see Ref. j^]). 
Because we do not use the definition of Ref. {4] , we explicitly check how our definition for the 
sign changes with N. For that purpose we fix /i/y/X ~ 0.95 and measure (cos#) for b = 10.00 
and N = 10, 20, 40. The results are presented on a logarithmic plot in Fig. |3] and confirm that 
the sign drops approximately exponentially with N. 

Before we proceed we wish to remark on the way the average sign behaves for \x < m^/2. 
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FIG. 3: A logarithmic plot of (cos(0)) versus N at b = 10.00 and fi/y/X ~ 0.95. 



Chiral perturbation theory gives different predictions for this behavior, depending on the value 
of the quark mass, the chemical potential and the volume. For example, in the thermodynamic 
limit, taken at fixed quark mass, the sign should freeze at one for all \x < m n /2 28j. In 
the so-called 'epsilon regime', where the quark mass and fi are decreased when one takes 
the thermodynamic limit, the average sign approaches a smooth function of the variable 
z = (2fi/m n ) 2 , that is equal to unity at z = and to zero at z — 1. Therefore, interpreting 
the large- N limit as the thermodynamical limit of our system, we expect that the large-iV 
limit of of (cos#) will be 0(1) for z < 1 and zero for z > 1. From Figs. ([2]H© it is fairly clear 
that for z > 1 our results are consistent with this expectation. The situation for fi < m n /2 is 
less clear since we do not know whether our quark mass and the values of [i that we study are 
inside the epsilon regime or not (for that we would need to measure /„-). Nonetheless, the fact 
that the sign clearly drops with iV suggests that we are either inside or close to the epsilon 
regime in that part of the phase diagram. Thus, as a guide to the eye we plot in Fig. |5] the 
average sign in the epsilon regime for our system. Here we assume that we are close enough to 
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the continuum limit, and so we plot the quenched averaged sign (e 4ie> ) quenc h e d- 8 The factor of 
4 in the exponent reflects the fact that we use naive fermions in our numerical studies; in the 
continuum limit these fermions are four-fold doubled which means that our Dirac determinant 
det (Aattice) factorizes to (det D Nf= i^j , where D Nf=1 is the Dirac operator of a single flavor 
theory in the continuum limit. 

VII. TESTING Z N AVERAGING. 

We now proceed to test the averaging procedure by calculating the average sign of 
(det.D) z for the same parameters presented in Fig. El We denote this quantity by (cos9)z, 
although we emphasize that what is measured is the sign of (det D) z and not the average 
of cos 9. 

We present the results in Fig. H] where we see that Z^- averaging makes a dramatic differ- 
ence: throughout the range of chemical potentials G [0, 1.2] the .Z^r- averaged quantity 
(det.D) z is real and positive. In fact, the N dependence of {cos9)z suggests that this pos- 
itivity will persist until about — 1.25, a regime that includes, as a subset, all of the 
hadronic phase. This is most clearly seen by comparing Fig. H] with Fig. [2] — there is a wide 
regime with /i > m n /2 where without Z^ averaging the sign is exponentially small but with 
Z N averaging it is equal to one. Thus we see that, in that range, Z N averaging solves the sign 
problem of our theory. 

The absence of sign fluctuations below [ij \/\ — 1.25 does not mean that there are no sign 
fluctuations for all /i. In fact, for SU(40) there is a severe sign problem beyond the hadronic 
phase, in the range n/\f\ G [1.8, 2. 4]. 9 This is a real sign problem (i.e. it is not 'just' a silver- 
blaze problem) because it appears when observables start to be /i-dependent. Also note that 
this sign problem precedes the lattice saturation of the EK theory, where the sign fluctuations 
goes away again. 

8 I.e. we set p = 2 in Eq. (48) of Ref. 0] 

9 Our statistical errors in Fig.|4]are larger compared to those presented in Fig.[2]since we were able to calculate 
(cos 0) z for a fraction of our configurations in the SU(A0) case (recall that for each configuration we need 
to evaluate det-D for each of the 40 2 = 1600 terms in the sum in Eq. (|4.1|Q . Thus, our resources allowed us 
to only do so for 5000 of the 105000 configurations. 
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Average sign, b=10.00 




0.5 1 1.5 2 2.5 

FIG. 4: The Z^-averaged sign, (cos(0)) z , for b = 10.00 and N = 10,20,40. 



What happens if we only perform a partial averaging (see discussion at the end of 
Section H*V|) ? We present a comparison of the two type of averaging for N = 40 and b = 10.00 
in Fig. [5J The fact that the two averaged signs are nearly the same (contrast their behavior 
with the non-Z/v-averaged sign in the figure) means that a partial averaging works almost 
as well as the full one, and indeed we shall use it when our resources cannot accomplish the 
latter. An interesting point appears when we compare the results of partial Zn averaging for 
various values of N in Fig. [HI There, we see that, similarly to the full Z^ averaging case, the 
sign grows towards 1 for /i/\A £ 1.25 and fi/y/X Z 2.75. This means that in these ranges of 
//, partial Zn averaging removes the sign fluctuations from the averaging of det D in a way 
which is as good as the one provided by full Z^ averaging. 

Let us now ask whether Z^ averaging actually improves the signal to noise ratio of physical 
observables. Since unquenched observables always involve the calculation of (det-D), we ask 
what is the improvement that Zn averaging has to offer in the calculation of the latter (or 
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SU(40), b=10.00 




FIG. 5: The average sign with different versions of averaging. Here b = 10.00 and the gauge group 
is St/ (40). 

more precisely in the calculation of the free energy F = jj log(det-D)). 

To answer this question we check which of the following methods gives the smallest statis- 
tical error on F, provided we fix the computational effort. 

(1) No Zn averaging. Here we average a set of M gauge configurations. 

(2) With partial Zn averaging. Here we use a set of M/N gauge configurations. Recall that 
for given gauge configuration, we calculate the determinant iV times (see Section HV|) . 

We perform our check for N = 40, M = 40000, and b = 10.0 - a case where the sign problem 
can be relatively severe. The results of this check are given in Fig. [7J In the upper panel 
we show F, as obtained with both methods, and for convenience we show the corresponding 
average signs in the lower panel of the figure. A quick look at the figure shows that Z^ 
averaging can be useful. In fact, note that we do not show the results for F from method 
no. (1) when n/y/X > 2.25. This is because the severe sign fluctuations that appear in that 
method make the average of det D negative. In contrast, the average sign of method no. (2) is 
away from zero in the same regime and so (det D) pavtial z is positive and real there. The most 
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FIG. 6: The average sign with partial Zj\t averaging. Here b = 10.00 and the gauge groups are 
£17(10), SU{20), £17(40), and SU(60). 

important point, however, is that for fi/y/X < 1.25 and n/y/X > 2.75, the sign is exponentially 
small in N with method no. (1) while it is increasing towards unity with method no. (2). Since 
the numerical cost of method no. (2) is linear in N, we see that it provides a gain which is 
exponential in A7 10 To show the effect of the errors in the averages of det D on unquenched 
observables like ('ipip), we plot the results obtained for (ipip) with the two methods in Fig. [HJ 
Clearly Zjyr averaging is useful here. 

In the next section we show more results for F and (ipip)- These were obtained with Z^ 
averaging and for much larger statistical samples than the ones we discuss above. 

10 The fact that in method no. (2) one has less independent configurations might mean that its statistical error 
can be larger, but at sufficiently large- N and a moderately large number of configurations, the sign problem 
will always make method no. (2) preferable. This is already seen for 57/(40) in Fig. [7] in the range between 
fi/y/Xe [1,1.25]. 
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FIG. 7: A test of the efficiency of averaging for SU{4Q) and b = 10.0. Squares(blue) show the 
result of method no. (2) while circles(red) those of method no. (1). In the upper panel we show F 
and in the lower one the average sign. 

VIII. THE FREE ENERGY AND THE CHIRAL CONDENSATE 

In this section we present several results of physical interest, beginning with the way the 
free energy F behaves as a function of \i. We calculated F using averaging and present 
the results for b = 6.0, 10.0 and iV = 10, 20 in Fig. |9j Performing Z^ averaging for N = 40 
was too challenging for our resources and so in that case we only used partial Z^ averaging. 
Also, while the S77(40) results are similar to those of 577(10) and SU (20) for /i/vA < 1.5, the 
errors in the regime of the 'true' sign problem make the data not useful there - this is why 
we do not present the data for b = 10.00 in the regime n/y/X G [1.5,2.75]. Comparing this 
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FIG. 8: A test of the efficiency of Zjy averaging for SU(40) and b = 10.00 for the chiral condensate. 
Squares(blue) show the result of method no. (2) while circles(red) those of method no. (1). The 
data obtained with method no. (1) is slightly shifted on the x-axis, to distinguish it from the data 
of method no. (1). 

figure to the upper panel of Fig. [TJ we see that, beyond the hadronic phase, when the true 
sign problem is at its peak, even Zn averaging does not help. This is reflected by the fact that 
while the smaller statistical sample used to produce Fig. [7] gave positive values of (det D), the 
statistical sample we use in this section, which is 100 times larger, resulted in negative values 
of (det D) in the same regime. 

There are two additional important points to take away from Fig. [91 First, we see that 
the free energy is completely independent of \x throughout the hadronic phase, as we expect 
physically for T = 0. This is another way of seeing that the silver-blaze problem is absent 
from our calculation. Second, we see that F becomes ~ 2/i for large value of \i. This is the 
lattice saturation we discussed above (see Eq. ( 16. 2ft ). 

To emphasize that the results in the preceding sections are truly unquenched we present the 
behavior of the quenched condensate ($1/1) quenched in Fig. [KB The sharp change in ($1/1) que nched 
at around m vr /2 is the same as seen in other quenched lattice calculations, and it is a non- 
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FIG. 9: The Z^-averaged free energy density F z (see Eq. <^M>) of SU (10), SU (20), and SU (40) 
at b = 6.0 and 10.0. The calculation for the SU(40) case was done with partial Z/v-averaging. To 
guide the eye we connect the data with solid lines. 

sensical result reflecting the mutilation that the quenched prescription causes to the gauge 
theory at nonzero fi. Comparing Fig. [10] to the plot of the free energy, Fig. [9] we see that the 
drop in the former happens at values of /i where the free energy is still independent of fi. 
We attempted to use Zn averaging and calculate the unquenched condensate, but for 
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FIG. 10: The quenched quark condensate for N = 10, 20, 40 and b = 6.00, 10.00. 

S77(40) this proved too costly for our resources. 11 Instead, we calculated it using partial Z N 
averaging (see Section HV|) . Let us now compare the quenched and unquenched condensates. 
We begin at /i = 0, where we expect both condensates to be equal when N = oo. We test this 



11 The bottle neck was not the generation of configurations, but rather the repeated calculation of the observ- 
ables for each gauge configuration. For example, analyzing 10 4 measurements of SU(A0) for a single value 
of fi with Zn averaging would take around 555 hours on a 2.66GHz CPU. 
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in Fig. ITT1 where we plot the difference (^(?/>?/>)quenched — (V'V'} unquenched) /N for the three lattice 
couplings 0.60, 6.0, 10.0 and versus 1/N. As the figure shows, while the difference indeed goes 
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FIG. 11: The difference between the quenched and unquenched chiral condensate for b 
0.60,6.00, 10.00 versus 1/N for y, = 0. 



to zero for b = 0.60, it does not seems to do so for b = 6.0 and 10.0. We interpret this as 
resulting from working at too small values of N. Similar differences between quenched and 
unquenched averages at /i = were also reported for the overlap operator in this system 291 ]; 
it would be useful to understand this slow convergence to N = oo, but this is beyond the 
goals of our current study. 

We now turn to compare the quenched and unquenched condensates at \i ^ 0. In this case, 



our discussion in Ref. 



301 ] shows that the quenched approximation fails and that we should 



expect large 0(1) deviations (that do not go to zero at large N). While we explain this failure 
in Ref. [30j let us think about it from the point of view of the numerical calculation. If the 
quenched approximation was exact at large- N it would mean that the operators tpifi and det D 



2s 



are classical operators that obey large- N factorization since then we would have 

x detD) N =°° (^) x (detD). (8.1) 

At nonzero /x we know that Eq. ( 18. 11) fails because quenching fails — see the fictitious phase 
transition at /x — m n /2 observed in Fig. ([HI]) . Thus we expect correlations generated by the 
strong fluctuations that both ifjifj and det D have within the complex plane. Put differently, 
the angles 9 and a defined through Eqs. ( I3.18143TT9|) now fully spread over the range [—it, tt) 
and can become correlated even when N = oo. We will investigate this issue numerically in 
the next section. 

A direct measure of such correlations is the way the difference between the quenched and 
unquenched chiral condensates behave as a function of /x. We already saw that at the values 
of iV that we work with, the /i = condensates are different. This, however, we associated to 
significant 1/N corrections and distinguishing these from the 0(1) differences that we expect 
at nonzero /x, is hard to do unambiguously. Nevertheless, we attempt to do so by presenting in 
Fig. [12] the way the quenched and unquenched condensates, normalized by their fi = value, 
change with /x, for b = 10.0 and iV = 20 and 40. 

There are a few useful observations we can make on Fig. [12j First, we see that the behaviors 
of the quenched and unquenched condensates as a function of /x strongly differ for chemical 
potentials that are above half the pion mass (denoted by the vertical magenta dashed line), as 
expected from chiral perturbation theory. Specifically, we see that despite the sharp change 
in (ijjip) quenched at /x ~ m n /2, the physical unquenched condensate is unchanged throughout 
the hadronic phase and beyond. Second, we see that, within our statistical errors, the results 
for N = 20 and 40 are nearly on top of each other. Thus, combined with the theoretical 



expectations of Ref. [30|], it does not seem likely that the differences between the quenched 
and unquenched averages that we see for /x > m n /2 are 1/N effects. Finally, we see that the 
errors on the unquenched condensate greatly increase in the regime of the real 'non-silver- 
blaze' sign problem (see Fig. |3]). All these facts are in agreement with physical expectations. 
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Quenched and non-Quenched chiral condensate, b=10.0 
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FIG. 12: The quenched and unquenched chiral condensate, normalized by their = value, versus 
/i. Here b = 10.00 and we present results for iV = 20, 40. The vertical magenta dashed line denotes 
our estimate for half the pion mass (see Table [I]) . 

IX. DISTRIBUTIONS IN THE COMPLEX PLANE 



This section has three goals. First we wish to examine the correlations between ifjifj and 
det D that cause the failure of the quenched approximation. Second, we wish to see if our 
results for the distribution of 9 (as defined by Eq. (I3.18P ) varies with \i as predicted by 
chiral perturbation theory Q]. In that reference, the authors showed that for /i < m n /2, 6 is 
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distributed as a periodic Gaussian, while for fi > m n /2 the distribution becomes a periodic 
Lorentzian. Also, the widths of these distributions go to infinity in the thermodynamic limit 
(which in our case is the large- N limit). Third, we wish to see if the distribution of the 
eigenvalues of D in the complex plane is consistent with chiral perturbation theory. 

We begin with analyzing the correlations between the phase a of ipip and 9, and present 
these for SU (20), b = 6.0, and different values of [l in Fig. H3J What we learn from the figure 
is that for most values of fi there is some correlation between a and 9. Especially, it seems 
that the fluctuations of a and 9 are not small, and that for \i > m^/2, they spread in the 
range [— tt, 7r). We note that this, however, does not mean that they are strongly correlated 
since their joint probability distribution can still be approximately separable. 

The probability distribution of 9 itself is shown in Fig. [141 where we see that, as Ref. jj] 
anticipates, it is a periodic Gaussian for low values of fi. At larger values of \i it grows wider 
and flattens. We fitted P(9) with the forms 

1000 

Pg(9) ~ E ex P 

n=-1000 

P L (9) ~ l/[A-cos0] (periodic Lorentzian). (9.2) 

We present the resulting parameters of the fits in Table [Til The statistical error on the fit 
parameters represents the one-sigma error from the fit, as well as the differences between fits 
done for histograms generated by cutting the full statistical sample into N hin = 15, 25 and 35 
bins. The different values of the x 2 i n t ne table represent this dependence on iVbin- Our fits 
were uncorrelated, and changing iV^n is the way we check for the correlation of the different 
bins in the histogram. 

From the table we see that while the periodic Gaussian can be fitted to the data for low 
and high values of /i, the Lorentzian form cannot. At ~ 1.7-y/A the Lorentzian becomes 
consistent with the data, but at that point the histogram is so flat that both the Gaussian 
and the Lorentzian can be good fits. Note also that the values of x 2 are n °t too small, but 
looking at the plots they seem to reflect a general scatter of the data around the curve and 
so are likely to come from an underestimate of the statistical error. 

Let us now ask how Fig. [13] and Fig, [ill] change when we increase N. The results for SU(A0) 
are given in Figs. [T5lfT6l and we see that increasing iV makes the distribution of both 9 and a 
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27m) 2 /A (periodic Gaussian) 



(9.1) 



fi/y/X Gauge group 


Gaussian fit 


Lorentzian fit 


A X 2 /dof 


A X 2 /dof 


SU(20) 

0.31579 


0.71(3) 1.6-2.7 
1.17(1) 0.8-1.6 


No fit 
No fit 


SU(20) 

0.78947 

SU(4Q) 


6.5(1) 1.8-2.5 
11.5(5) 0.9-1.3 


No fit 
9.0(7) 0.7-1.0 


SU(20) 

1.7368 

SU(4Q) 


17(2) 1.8-3.1 
18(3) 1.2-1.5 


40(20) 1.8- 3.0 
50(30) 1.6 - 1.5 


SU(20) 

2.8421 

SU(40) 


2.84(4) 1.0 - 2.8 
4.65(5) 0.6-0.9 


No fit 
No fit 



TABLE II: Results of fits to the histograms in 6 (see text). 



more spread out. The corresponding fitting results to the distribution P(9) in this case are 
also presented in the Table [TT1 

It is hard to conclude from the figures how the correlations between 9 and a changes with 
N. Nevertheless, the widening of the histograms in Fig. [16]compared to Fig. [T3]is qualitatively 
consistent with Ref. |4| if we interpret the large- N limit as the thermodynamical limit. More 
accurate studies are required to pinpoint the nature of the correlation between ipip and det D 
(the absolute values of ipip and detD may play a role also). Probably the best indicator for 
this is the connected correlation function between these operators, which is nothing but the 
difference between the unquenched (ipip) and (V^quenched- As we mentioned in Section |VIII[ 
this difference is significantly nonzero and independent of iV for fi > m^/2. 

We note in passing that we see that the x 2 f° r S 11(40) cases are lower than for SU (20) 
and this may reflect the presence of significant 1/N corrections that are required to make the 
the fitting ansatze in Eqs. I9.1H9.2I consistent with our data. Also, we see that, as predicted 
by chiral perturbation theory, the Lorentzian is consistent with our data at /i ~ 0.8^X. This 
should only be seen as partial evidence to support these predictions since for the same values of 
/i, the Gaussian is also consistent with our data, and because for SU(20) it is only the Gaussian 
that provides an acceptable fit at that value of fi. Finally we note the way A depends on is 
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FIG. 13: The correlations between a and (see text) for SU (20), b = 6.0, and various values of \x. 
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FIG. 14: The probability density of 6 for 5C/(20), & = 6.0 and various values of \i. The solid curves 
are the Gaussian fits (see Table Hi]) . 
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FIG. 15: Same as Fig. US but for SU(40). 
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FIG. 16: The same as Fig. [T^l but for SU(40). The only panel that includes both the Gaussian 
and the Lorentzian fit is that of [i ~ 0.79-^/A. The Gaussian curve is a solid line (blue), and the 
Lorentzian is the dashed (red). It is clear that both are very close to each other. 

expected to be linear for the Gaussian case (see Ref. but this is only partially consistent 
with what we see in Table [TTJ this is expected to happen only for /j, ~ 0.32y/X < m 7r /2, but 
there we see that Asu(40)/^su (20) = 1-65(7) instead of 2. In the case of fi ~ 2.8^, where 
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P(0) is clearly Gaussian and not a Lorentzian, we see that A S u(4o)/ '^su{20) — 1-64(3). We 
do not fully understand these results which actually correspond to Asu(N) ~ jy - 71 ® -0 - 72 ^. 
They may reflect significant 1/N corrections in the SU (20) case. 

We end this section by presenting the way the eigenvalues of D are scattered in the complex 
plane. According to chiral perturbation theory, when fi < m n /2, the bare quark mass m is 
outside the support of the eigenvalue distribution, while for /x > m^/2 it is within it. This 
is confirmed in Fig. [T7J where we show the eigenvalue scatter of 50 gauge configurations for 
SU(40) and b = 6.00. We present the cases of fi/y/X = 0.31579, 0.60, 0.78947 which are below, 
close to, and above \ra^/y\. Also, the cases of fi/y/X ~ 0.31579, 0.78947 correspond to the 
two upper panels of Figs. [ToT - TlBl Recall that in our current work we fix the quark mass to be 
m/VA = 1/Vtt ~ 0.5642 (see Eq. (ETJ). 



X. SUMMARY AND CONCLUSIONS 



In this paper we studied the sign problem of the euclidean path integral of large- A QCD 
at nonzero chemical potential \i. To do so we used the Eguchi-Kawai (EK) equivalence which 
allows us to approach the problem via the volume reduced version of theory. For reasons of 
computational cost we focused on the two dimensional case in this paper. Extensions to four 
dimensions are straight forward and we discuss them below. 

We used lattice Monte-Carlo simulations to explore the sign problem numerically and did 
so for three lattice spacing, and for SU(N) gauge groups with A < 60. We measured the 
average sign of det D, and saw that it decreases to zero when // ~ m n /2, but then rises back 
to 1 when the lattice site is saturated with quarks. 

We proposed a way to suppress the sign fluctuations in the hadronic phase which amounts to 
replacing the pure-gauge average of det D, by an average over a real and positive quantity. We 
denote the latter by (det D) z , and it is the average of det D over a set of gauge configurations 
that are related by center transformations. We call this method '^TV-averaging' (see 
also [8]) and we test it numerically. We find it to be successful and that, in the large- A 
volume-reduced system, the average sign of (detD) z is 1 for a wide range of \i that includes 
all the hadronic phase. Thus, this method removes the sign fluctuations from the hadronic 
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FIG. 17: The scatter in the complex plane of the eigenvalues of D. Here b 
group was SU (40). 



6.00 and the gauge 



phase. The computational cost of averaging grows linearly in N, and this means it 
provides an exponential gain in calculating the average of detD, since the computational 
cost of overcoming the sign fluctuations involved in such an average by brute force grows 
exponentially in N. This means that averaging solves the silver blaze problem in our 
volume-reduced large- N theory. 

Z^-averaging also identifies the regime in \i where the true, non-silver-blaze, sign problem 
of the large- N reduced theory is most severe. The latter happens beyond the hadronic phase, 
just before the saturation regime. We identified the saturation regime also in measurements 
of the free energy, which becomes linear in \x in that regime. Our measurements of the free 
energy and of the chiral condensate are unquenched and as anticipated physically, we see 
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that they show //-independence throughout the hadronic phase. Our measurements of the 
quenched chiral condensate serve to contrast this since they show that (ipip) quenched changes at 
around /i ~ m n /2, where nothing happens to the free energy and to the physical unquenched 
condensate. This reflects the breakdown of the quenched approximation. 

In an attempt to numerically understand the failure of the quenched approximation men- 
tioned above, we analyzed the way the operators ipip and det D are correlated with each other 
for different values /i. Specifically, we checked whether their correlations is mostly due to a 
correlation between the phases of the operators, but were not able to unambiguously conclude 
if this is the case or not. Instead, we directly calculated the connected correlation of ipip and 
det D. This correlation is nothing but the difference between the quenched and unquenched 
chiral condensates. We saw that this difference becomes large for fi > m n /2, and that it does 
not go down with increasing N. This means that the failure of the quenched prescription can 
be thought of as the breakdown of large- N factorization. It would be useful to understand this 
breakdown from a physical point of view rather than a numerical one. It will also be useful to 
understand our results for the difference between the quenched and unquenched condensates 
at fi = and verify that they go away at iV — > oo by performing simulations for large values 
of N. 

We also measured the way the phase of det D is distributed as a function of iV and /i, and 
saw that, in agreement with chiral perturbation theory (see Ref. [4]), these distributions are 
Gaussian for fi < m 7r /2 and become extended when \i increases beyond that. The width of 
these distributions increases towards the thermodynamical limit of iV — > oo. At the quantita- 
tive level, our results support the predictions of Ref. {4] only partially, but this may be due to 
1/N effects. Finally, we confirmed that when > the quark mass m enters the support 

of the density of eigenvalues of D. 

The analysis we presented here can in principle be repeated for four dimensions but this 
will require more work. In particular, in 4c?, the straight forward EK theory is not equivalent 
to large-iV QCD (even in the hadronic phase) and other prescriptions, like those suggested 
in Refs 12] and 13] (and studied numerically in Refs. [yj] and 15]), are needed to overcome 
this issue. Once this is done, and one has a set of pure-gauge field configurations in hand, 
then the methods of this paper can be used with a modest increase in computational effort. 
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Specifically, the only change is that the dimension of the Dirac matrix D is doubled relative 
to the two-dimensional case. We therefore estimate that future 4c? studies of the issues we 
explore in this work are feasible. 

Finally, we proposed that averaging solves the silver blaze problem when the contri- 
bution to the Dirac operator determinant from baryonic worldlines (or Polyakov loops whose 
winding number is a multiple of N) is suppressed on a configuration-by-configuration basis 
(for an explanation of this see Section IVBI) . In the current paper we see numerically that in 
1 + 1 dimensions such suppression indeed takes place in the large- N volume-reduced version of 
the QCD, but we do not know if this is also the case for physical (four dimensional and 3-color) 
QCD. More precisely, in Section IVBI we discussed how such a suppression will not happen 
at T — > since, in that limit, a configuration with k baryonic worldlines will most likely be 
only moderately suppressed by the Boltzmann factor of e - kNm ^/( 2T ) . We argued, however, 
that there may be a region of T > where, instead, these worldlines might be suppressed and 
where averaging would work. This 'golden region' is the analog of the 'golden window' 
in the studies of baryon correlation functions , where one searches for optimized baryonic 
wave functions whose noise couples very little to pion systems. Therefore it seems interesting 
to attempt Z^ averaging in physical QCD (with light quarks) and try and locate this golden 
region. 
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